Exact real-time dynamics of the quantum Rabi model 
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We use the analytical solution of the quantum Rabi model to obtain absolutely convergent series 
expressions of the exact eigenstates and their scalar products with Fock states. This enables us to 
calculate the numerically exact time evolution of {ax{t)} and {(Tzit)) for all regimes of the coupling 
strength, without truncation of the Hilbert space. We find a qualitatively different behavior of both 
observables which can be related to their representations in the invariant parity subspaces. 
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I. INTRODUCTION 

The quantum Rabi model (QRM) describes one of 
the simplest strongly coupled quantum systems: a single 
bosonic mode interacts with a two-level system. It forms 
the basic building block of theoretical approaches to the 
interaction of light with matter [1] or to electron-phonon 
interaction [2] . Due to recent progress in nanofabrication, 
the strong-coupling regime of the QRM could be accessed 
experimentally within circuit QED [3, 4], where the two- 
level system corresponds to a single qubit. To construct 
controllable solid-state quantum gates, it is necessary to 
comprehend the behavior of the QRM for all regimes of 
the coupling strength [5]. 

Although the lower part of the spectrum of the QRM 
can be computed numerically to arbitrary precision us- 
ing exact diagonalization on a truncated state space due 
to the convergence properties of the associated contin- 
ued fraction [6-9], it is by no means evident whether the 
properties of the numerically obtained eigenstates corre- 
spond to those of the true eigenstates of the untruncated 
model. Even if the spectra of the truncated and the full 
model are close for a subset of the eigenvalues, the corre- 
sponding eigenvectors could still be related by a unitary 
transformation, which does not necessarily approach the 
identity. Consequently, a variety of systematic approxi- 
mations making use of different basis sets in the infinite- 
dimensional Hilbert space have been developed [2, 10-15]. 
However, all of them implement a truncation procedure 
at some stage of the calculation and thus cannot answer 
the basic question on the relation of the truncated to the 
full model. 

Beside this theoretical goal, there is a quite practical 
reason to investigate the true eigenbasis of the QRM: To 
compute the dynamics of the system on all time scales, 
not only the exact eigenvalues, but also the exact eigen- 
states are required. Only then one may hope to pre- 
dict reliably the behavior of qubits within a circuit QED 
setup [5]. 

In this paper we solve the problem by means of the 
recently obtained exact solution of the QRM [16]. We 
construct absolutely convergent series expansions of the 



eigenstates in terms of known basis sets, which allows to 
estimate the true error in any numerical evaluation. It 
turns out that a standard calculation with double pre- 
cision, sufficient to compute the spectrum, fails for the 
eigenstates. 

The QRM Hamiltonian reads 
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The Hilbert space is given by L^(R)(8)C^, where the basis 
of the spin part of the wavefunction | cr) is {| -I- 1), | — 1)}. 
We employ the Bargmann space B of analytical functions 
to represent elements of L^(M) [17]. In B the creation 
(annihilation) operators a^ (a) are represented as z [dz) 
and normalization is defined with respect to the scalar 
product 
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An arbitrary analytical function (j){z) is an element of the 
Bargmann space B if {(f>\(j)) < oo. 

Each eigenstate of i^R belongs to one invariant sub- 
space (parity chain) [18] with a fixed Z2 quantum num- 
ber, the parity, taking values ±1. The parity operator 
for Hb, is defined as P == (-I)^^-(t^. Due to [TJr, P] = 0, 
the Hilbert space is a direct sum of invariant subspaces 
H+ ffi H- and H^ takes the form H± on 7i±. 

The invariant parity chains 'H± are spanned by 
{(t>s{z) (g) I ± 1), (/)a(2;) (^ I T 1)}, where (t)a.,s = (l/2)[0(z) T 
(j){—z)] is the (anti-)symmetric part of (j){z). Each parity 
chain is therefore isomorphic to B via the mapping 
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where (j)s{z)(E)\±l) + (f)a{z)(E)\^l) = 10, ±) is an element 
of 'H±. Using this isomorphism, the Hamiltonian H± 
reads 



H± = ujzdz + g{z + dz) ± AT 



(4) 



with A — wo/2, and the reflection operator T acts on el- 
ements of S as T{(p){z) — 4){-~z). For notational simplic- 
ity, we restrict ourselves to ■?{+ and measure all energies 
in units of lu. 



The eigenfunctions ipmiz) and eigenvalues Em of if+ 
have been obtained in Ref. 16 as two equivalent repre- 
sentations, 
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'ipmiz) =(/)2(x„,-z) ^ e^"" ^ Kn{x,n){-z + g)" , (5a) 
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where Xm — -E-m + 9^ and 

nKn{x) ^ .fn-i{x)Kn-i(x) - K^-iix), (6a) 

K^{x)^Ux), Ko{x)^l, (6b) 
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Both (5a) and (5b) converge for all z in the complex plane 
if and only if Xm satisfies the spectral condition [16] 



G+(xm,0) ~ where 

G+{x, z) = (l)2{x, -z) - (t)i{x, z). 
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For all other values of the spectral parameter x ^ x^, the 
functions 0i^2(a;, z) are defined via (5a, 5b) only within 
the circle 1^^ + g| < 23. Furthermore the eigenfunctions 
are not normalized when expressed in terms of (5a, 5b). 
Both difficulties are overcome in the following. 



II. TIME EVOLUTION 

If an initial state \(f)Q,+) with fixed parity p = 1 is 
prepared at time i = 0, the time evolution takes place 
within 'H+ and reads 



|0(i))=e-^^+*|</>o)=^e-^^i*|V„ 
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The task consists in computing (fAml'/'o) and with that 
the norms (V'mlV'm) for m = 0, 1, 2, . . .. 



A. Formal solution 

Starting from (5a, 5b) we obtain an expansion of |V'm) 
into normalized cigcnstates of the shifted harmonic oscil- 
lator, 
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as follows: 
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From these equations and the fact that the eigenbasis 
of the shifted harmonic oscillator {|n; g)} is orthonormal 
we obtain two equivalent representations of the squared 
norm of the eigenstates: 
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Using the expansions of the initial state 0o into shifted 
oscillator states 0o = Yln=o ^^\'^'^ ^ d) = Hn=o P^\'^'^ 9) 
we have 
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^V^i^„(x„)(-l)"a„, (13a) 
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The preceding equations provide a formal solution of the 
time evolution problem of (9). 



B. Asymptotic numerics 

In practice a numerical calculation of (12,13) fails as 
the sums diverge for all x which do not exactly coincide 
with a value Xm of the spectrum. For every approxima- 
tion Xm to Xm, the absolute values of the summands of 
the series in (12,13) start growing for an n > n', with n' 
unknown a priori. In Fig. 1(a), the dashed lines show this 
behavior for the example of the summands of the norm 
calculated via (12a). The picture suggests the numer- 
ical finding we confirmed for different parameters, that 
n' grows for growing eigenenergies Xm- We further find 
numerically that n' is large if the error of the approxi- 
mation 5 = \xm — Xm\ is small, which is characteristic 
for series that are asymptotically valid. It will be shown 
in the next section that this asymptotic validity holds 
rigorously true for the series in (12,13). 

Due to its asymptotic validity for (5 ^ 0, we can eval- 
uate (12a) by approximating 



JV 

E^ 

n=0 



(V'm|V'™) = e9 ^n\Kl{5:m) + R 



N- 



(14) 



Here the cutoff N must be chosen large enough so that 
the relative error Rn / (ipmlipm) , estimated as 
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is as small as desired. If the evaluation of (15) is done 
with a low-precision approximation x^ of x^, the se- 
ries might start to grow before the precision of 6^01 is 
achieved, i.e., n' < N. In this case the calculation has to 
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FIG. 1. Panel (a) shows the square root of the absolute values 
of the summands y/n\\Kn{xm)\ of (12a,14). Panel (b) shows 
the absolute value of the coefficients |7i'n(a;m)| (see Eq. (6)) 
of the eigenfunctions and panel (c) shows the relative error 
erGi(A'^ ~ n) (15) of the series expression in Eq. (14). All 
quantities are depicted for three eigenenergies {a;o, a^s, a;2o} 
and the parameters g/cj — 0.7, A/lj = 0.25. We present a 
calculation with double precision and 5 = 10"^" (dashed lines) 
and a calculation with high (50 digit) precision and S — 10^'^" 
(solid lines). The horizontal line in panel (c) marks the value 



be repeated with a better approximation Xm, i.e., with a 
smaller error S. The expressions (12b) and (13) can be 
treated analogously. 

In Fig. 1 we illustrate the asymptotic behavior of the 
sum in (12a). We want to evaluate it up to double preci- 
sion (i.e., by demanding e^ei = 10~^^) which is in general 
the required precision for possible further computations. 
We evaluate (12a) for two different values of Xm ■ The first 
value is calculated with double precision and 5 = 10^^°. 
Although a computation of the norm is possible for the 
lower part of the spectrum, for the twentieth eigenvalue 
X20 the error ej-oi does not become as small as the de- 
sired accuracy 10~^^, but starts to grow before. In a 
second step we evaluate the norm with an approxima- 
tion Xm obtained in a high (50 digit) precision calculation 
{S = 10~^°) [19]. This setup is depicted by the solid lines 
in Fig. 1(c) which cross the horizontal dashed line, mark- 
ing Crci — 10~^^, also for 2:20. However, the coefficients 
Kn{xm) in panel (b) are almost identical for the high- 
and low-precision calculations for n < n' ^ i.e., before eroi 
starts to grow. This observation serves as an a posteri- 
ori confirmation of the asymptotic validity of (12a). In 
the regime where the series is asymptotically meaningful 
\K„\ < ^= holds, while for n > n', \Kn\ decays expo- 
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of vn!. These regimes are separated by clear kinks in 
Fig. 1(b). 

The problems that arise when choosing a low-precision 
calculation are related to the structure of the three-term 
recursion (6) that determines the sequence Kn{xm)- For 
high values of Xm, the absolute value of Kn{xm) becomes 
so high that the sign-changing behavior of the sequence 
leads to large numerical errors. For this reason, in a 
double-precision calculation even the spectrum cannot be 
determined for a; > 40 as the computation of the function 
in Eq. (8) becomes numerically unstable. The results of 
Sec. IV are all calculated using the high-precision setup 
of Fig. 1. 



III. ABSOLUTELY CONVERGENT 
REPRESENTATION OF 4>i,2(x,z) 

The reason for the divergence of (12,13) is that 
4>ia{x, z) in (5) is ill-defined for |z 4- (?| > 2g if a; ^ Xm- 
In the following, we present an analytical continuation of 
01,2 (a;, z) to the left complex half-plane ^{z) < that is 
well defined for all x. Using this analytical continuation, 
we perform the integration of the scalar product (2) over 
the full complex plane and obtain finite values for the 
norm (12) and scalar products (13) also for x ^ Xm- For 
X = Xm this evaluation yields the values associated to 
the true eigenstates. A comparison with the results of 
the preceding subsection then establishes the asymptotic 
validity of formulae (12,13). 

The squared norm (i/'mlV'ni) reads 
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'''lp{Xm,z)')p{Xm,z). (16) 



To compute (16), one needs ^{xm,z) in the whole 
complex plane. Using the analytical continuation of 
01 (-z), 02 (2) into the left complex half-plane, we define 
for arbitrary real cc. 
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nentially, which is too slow to compensate the growth 



The function '0(x, z), defined by (j)i{x,z) in the left and 
by (f>2{x,~z) in the right half-plane, exhibits a disconti- 
nuity along the imaginary axis ii x — g^ ^ speci?+, but 
becomes analytic everywhere if a; = Xm- By replacing 
(16) by (17), we obtain finite expressions for all x and 
the correct ones ii x — g'^ belongs to the spectrum. 

The analytic continuation of 4>i^2ix,z) from the circle 
with radius 2g around z — —g to the left half-plane is 
constructed with a Mobius transformation. The system 
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FIG. 2. Image of left z half-plane (light blue) in the w plane 
using the Mobius transformation (19). 



of ordinary differential equations satisfied by 0i (z), 02(-2) 
reads [16] 

(z + g)<i,\ (z) + (5z - X + g')4>x{z) + A^sC^^) = 0, (18a) 
(z - g)0^(z) - (5z + X - 5')02(^) + A0i(z) = 0. (18b) 

Primes denote differentiation with respect to z. The sys- 
tem is characterized by two regular singular points at 
z = ±3 and one irregular singular point at infinity [20]. 
Consider now the transformation, z — )■ w(z), 

w{z) = — -, z{w)^g -, (19) 



z/g-S 



w — \ 



which maps the irregular singular point at z = oo to 
w = \ and the two regular singular points at z — ^g to 
w = and w = —1 respectively (see Fig. 2). The left z 
half-plane is mapped onto the disk with radius 2/3 and 
center w = 1/3 whereas the imaginary z-axis maps onto 
the boundary of the disk. With x — 4g^ — x, the system 
(18) reads in the variable w, 

wiw ~ l)^'i(«;) - (^ -I- x)4>iH = A4>2H, (20a) 

^^2M + (^ + 65' - i)U^) = A0i(u;). (20b) 

It has regular singular points at w = — 1 , and one irreg- 
ular singular point at w = 1, while w = (X) is a regular 
point. Now we may expand the functions (t>j(w) with 
J = 1 , 2 in a power series around the regular singular 
point w ~ 0, 
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(f>j{x,w) = Va(f)(a;)u; 
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and obtain, with the definition 
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for n > — 1 the following four-term matrix recurrence 

Vn+l = AnVn + B^Vn-l + C„t!„_2. (23) 



The coefficient matrices are given by 

/ {4g^ -2x+n){n-x)+2A^ 2A /-^ a+2n-2 n 
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The condition (f)j(x,z) = (j)j{x,w{z)) fixes the initial val- 
ues of Vq and Vi , namely 
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Together with v^i = 0, all t)„ for n > 2 follow recur- 
sively. 

The asymptotic behavior of the a„ can be inferred 
from (23) as lim„^oo |a„_[]^/a„ | — 1. Hence the radius 
of absolute convergence of the expansion (21) is 1 for all 
X, as expected from the singularity structure of (20): The 
expansion around w — converges up to the neighboring 
singular points, w — ±1. Now we can perform the in- 
tegral in (17) because the image of the left half-plane is 
contained entirely within the region of absolute conver- 
gence of the series expansion (21), which furnishes the 
desired analytical continuation of the functions (t>j{x, z). 
The integral (17) reads in w coordinates 

(H^mx)) = / ^^ V ^_.(.).M 

J n |w-l|^ 

disk 

x(\Mx,w)\^ + \M^,w)\^y (26) 

The numerical integration in (26) is done in polar coor- 
dinates over the disk in the w plane obtained by mapping 
the left half-plane with the Mobius transformation (19) 
(see Fig. 2). Fig. 3 shows that the integrand in (26) is 
a smooth function, for which the numerical integration 
is easy to control. The singularity that is still present in 
the wave function (j)j (z, w) for w = 1 is suppressed in the 
integrand by the regularizing factor e^^^'^'>^^'^\ As ex- 
pected from the convergence properties of (j)j{x,w), the 
representation of the norm in Eq. (26) depends smoothly 
on the argument x. This means that for values of x that 
differ only slightly from the spectral values used in Fig. 3, 
the corresponding plots of the integrand are essentially 
the same. 

Finally, we confirmed the validity of the asymptotic 
formula (14) by comparing it to a numerical integration 
of (17). The validity of the remaining expressions (12,13) 
was also checked. 
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FIG. 4. (CTz(f)) for A/o; = 0.25 and (a) ff/o; = 0.1, (b) gjuj = 
0.7, (c) gjbj = 2. The evolution starts with coherent states 
with a = (solid line) and a = 2 (dashed line). 



FIG. 3. The integrand of (26) for the first (a) and the fifth 
(b) eigenenergy (xq, = 0.06038... and xs = 4.9355...) in the 
positive parity subspace for parameters g/oj = 0.7 and A/o; = 
0.25. 



IV. SPIN OBSERVABLES 

In this section we calculate the time evolution of (p:^ (t)) 
and {(Jxit)) using the method presented in Sec. II. The 
"population inversion" ^^(i) has been thoroughly inves- 
tigated for the resonant case {g/uj ^ 1 and A/a; ~ 0.5) 
within the rotating-wave approximation, where it shows 
the typical collapse and revival oscillations [21] as well 
as within analytical [14] and numerical [22] treatments 
of the full model. Here we follow a different path by 
investigating the problem for the non-resonant case and 
arbitrary parameters in all coupling regimes. We find 
a qualitatively different behavior for {az{t)) and {(Jx{t)) 
for strong coupling. As will be discussed below, this dif- 
ference can be understood in terms of the Z2 symmetry, 
manifested in the two invariant parity chains. 



In terms of parity chains it reads, 

|^(0))=e-"'/2(|ch(az),+) + |sh(az),-)). (28) 
Using the isomorphism (3), we find for the time evolution 

of (<T.) 

{a,{t)) = e-"'[(cosh(az)|e-'^+'fe*-f^+*|cosh(az)) 

- (sinh(az)|e-*-f^-*fe*^"*|sinh(az))], (29) 

where T is the reflection operator in B. This expecta- 
tion value is plotted in Fig. 4. As can be seen from (29) 
the time evolution is solely governed by the commutator 
[T, H±] and the system evolves separately in each par- 
ity subspace. Thus, when evaluating the expression, one 
can assign a definite parity to the frequencies Ef^ — E^. 
As these frequencies are always of order (n — m)uj for 
sufficiently strong coupling (7/A, the behavior on time 
scales of a few 2tt/uj is dominated by oj as the smallest 
frequency. 



A. Time evolution of (az) 

To study the time evolution of a^, we prepare the 
system at time t = in a product of a coherent state 
I a) = e~" /2-i-Qz ^ii]j g^yj eigenstate of CTz 



|V(0)) = |a)®| + l). 



(27) 



B. Time evolution of (ax) 

In complete analogy to the case of <Tz , we prepare the 
system at time t = in a product of a coherent state 
with an eigenstate of ax 

1^(0)) = |a)®-^(| + l) + 1-1)). (30) 
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FIG. 5. {ax{t)) for the same setup as in Fig. 4. Note that (c) 
shows that {ax{t)} ~ 1 for aU times if the initial state is given 
by the hnear combination of two cat states in Eq. (32). We 
observe a sharp step-like decay (solid line) in (c) for a system 
prepared in the vacuum (a = 0). For a system prepared in a 
coherent state with a = g/u = 2, & constant plateau is seen 
in (c) (dashed line). 



Using the representation in the two parity subspaces, we 
obtain the time evolution of (cx) as 



^-iH+tiH-t 



{(^x{t)) = -(a|e --^e^ 



e-'"-'e'"+'\a) 



(31) 



This expectation value is plotted in Fig. 5. When eval- 
uating the expression in (31), the frequencies have the 
form E~ — E^ with eigenenergies from different parity 
chains. These frequencies may become arbitrary small. 
In particular the frequency E^ — E~ is small for strong 
coupling, because E^ and E~ both approach the same 
value nuj — g^ jio for growing 5/A from above and be- 
low. Therefore small frequencies dominate {crx{t)) for 
sufficiently strong coupling. We note that the frequency 
_E+ — E^ vanishes if the parameters g and A satisfy the 
Juddian condition Kn{nuj) = 0. 

For g/uj = 2 and A/uj — 0.25 (i.e., in the so-called deep 
strong coupling regime [18]) we observe two features in 
Fig. 5(c). The first is the sharp step-like decay of {ax{t)) 
for a system prepared in the vacuum (a = 0) at t = 0. 
Each step has the length 2tt/uj. This structure disappears 
for higher values of a. 

Secondly, {oxit)) stays perfectly constant for a — gjio 
as depicted by the dashed line in Fig. 5(c). In this case, 
the representation of the initial state (30) in the parity 
subspaces is a superposition of the two "Schrodinger cat 



states" |C±) e n± [23] 



|v[/(0)) = ^(|C+) + |C-)) for a = ^. 



(32) 



\C- 



.-9V2<^ 



(ch(gz/a;) (g) | + 1) + sh(.gz/a;) «) | - 1)), 
|C_> =e-^'l'^^' (ch(gz/cj) ® I - 1) + sh(gz/cj) ® | + 1)) . 

The states \C±) are eigenstates of i/± for A = and 
one would expect their approximate conservation over 
time for small A but decay at longer time scales or ap- 
preciable ratio A/(7. It is therefore surprising that the 
dynamical protection of \C±) (and concomitantly a^) ex- 
tends to long times (see next section) and applies to quite 
large A/g, which is 0.125 in Figs. 5(c) and 7(c). This 
remarkable feature of the QRM at couplings gjiij > 1 is 
related to the excellent approximation of the true ground 
state in each parity chain by the shifted vacuum state (see 
Eq. (10)) and will be discussed in more detail elsewhere. 
The dynamical protection of |^(0)) leads to the constant 
plateau (dashed line) in Fig. 5(c). Fig. 7(c) shows that 
this behavior persists even on long time scales. 



C. Long-time behavior 

In the preceding subsections we pointed out that for 
sufficiently high g/A the frequencies appearing in the 
time evolution of {az(t)) have an approximate lower 
bound of w while this is not the case for {ax{t)). This is 
inferred from an analysis of the representation in the par- 
ity subspaces and can be seen by comparing Figs. 4(b), (c) 
with Figs. 5(b),(c). Only in the latter slow oscillations 
are visible. 

This effect is also observed in the long-time behavior 
displayed in Figs. 6 and 7. For weak coupling {g/uj = 0.1) 
the time evolution of {az{t)) and {ax{t)) is rather similar 
in Figs. 6(a) and 7(a). For strong coupling {g/uj > 0.7) 
the same analysis as for short time scales applies. In the 
case of CTz , uj remains the smallest frequency and the os- 
cillations are dominated by higher harmonics, although a 
modulation by small beat frequencies, due to the super- 
position of contributions from both parity subspaces (29) , 
is visible. By contrast, the overall behavior of {ax{t)) 
for g/io > 0.7 is dominated by oscillations with period 
T* > 2tt/uj, as seen in Fig. 7(b) and (c). According to 
the argument of the previous subsection, T* grows with 
g. For g/uj = 0.7 the dominant frequency 1/T* is then 
essentially the same for initial coherent states with a = 
and a — 2, see Fig. 7(b). This is true also for the case 
of g/u! — 2, where it is seen that the step- like struc- 
tures follow an overall oscillation with period T*. Only 
in the special case of a = g/ui, in which the time evo- 
lution is dynamically protected, the initial expectation 
value is maintained for arbitrary long times (dashed line 
in Fig. 7(c)). 
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FIG. 6. {o"z(t)) for times up to 40 fundamental periods (To — 
2tt/uu) and the same parameters as in Fig. 4. The case g/uu — 2 
is spht into panels (c), a — 0, and (d), a = 2, for clarification. 
The effect of the modulating beat frequencies is not visible for 
g/uj = 0.1 on this time scale, whereas it appears clearly for 
g/uj — 0.7 and 2. In all three cases higher harmonics dominate 
the time evolution for a — 2 (dashed lines in (a) and (b), resp. 
panel(d)). 



V. CONCLUSION 

We have described a method to compute exactly all 
eigenstates of the quantum Rabi model together with 
their norms and scalar products with states of the stan- 
dard basis. In this way the development of an arbitrarily 
prepared initial state can be obtained for short and long 
times. The method uses the Z2 invariance of the model, 
the representation of the oscillator degree of freedom in 
the Bargmann space of analytical functions and a con- 
formal mapping to obtain absolutely convergent expan- 
sions of the norms and overlaps with Fock states. The 
method overcomes the limitations of numerical proce- 
dures or analytical approximations based on a truncated 
Hilbert space. Moreover, these techniques themselves can 
now be mathematically justified by comparison with the 
analytical solution and retain in this way their usefulness 



for practical calculations. This will be discussed else- 
where. The splitting of the full space into two invariant 
subspaces results in a qualitatively different behavior of 
(cr^(t)) and {ax(t)) on intermediate and long time scales: 
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FIG. 7. {o-x{t)) for times between and 40To as in Fig. 6 and 
parameters as in Fig. 5. The qualitative behavior depends 
clearly on the initial state for small coupling g/oj = 0.1 (top 
panel) but not for g/u) = 0.7 (b). The dominating period 
T* > To grows with g due to decreasing differences E^ — E^ . 
For a = g/co — 2, (ax) stays constant. 



{az{t)) is always dominated by integer multiples of the 
photon frequency ui. These oscillations are modulated by 
small beat frequencies because the QRM spectrum is not 
equidistant and both parity chains interfere. {ax{t)), on 
the other hand, shows oscillations with a dominant fre- 
quency much smaller than w, which is directly related to 
the mixing of both parity chains effected by the operator 
ax- The oscillation period grows with growing coupling 
between qubit and radiation field. For initial coherent 
states fulfilling the condition (32), we observe "dynami- 
cal protection" in the deep strong coupling regime which 
could be relevant to the control of fast quantum gates 
within circuit QED. 
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